%% process_showbasis
function process_showbasis(basis_type, derivative_type, varargin)
[x, y, z] = meshgrid(-1:0.05:1, -1:0.05:1, -1:0.05:1);
V = [x(:), y(:), z(:)];
phi = basis_ref(V(:,1), V(:,2), V(:,3), basis_type, derivative_type);
figure("WindowStyle", "docked");
hold on;
for i = 1:size(phi,2)
    idx = phi(:,i) == 1;
    scatter3(V(idx,1), V(idx,2), V(idx,3), 10, phi(idx,i), "filled");
end
hold off;
xlabel("X");
ylabel("Y");
zlabel("Z");
axis("equal", "tight");
view(120, 30);
end

%% basis_ref
function N = basis_ref(x, y, z, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 10);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y-z;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = z;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = 0;
    elseif strcmp(derivative_type, "dz")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 0;
        N(:,4) = 1;
    end
elseif strcmp(basis_type, "P2")
    lambda = 1-x-y-z;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda-1);
        N(:,2) = x.*(2*x-1);
        N(:,3) = y.*(2*y-1);
        N(:,4) = z.*(2*z-1);
        N(:,5) = 4*y.*z;
        N(:,6) = 4*x.*y;
        N(:,7) = 4*x.*z;
        N(:,8) = 4*y.*lambda;
        N(:,9) = 4*z.*lambda;
        N(:,10) = 4*x.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1-4*lambda;
        N(:,2) = 4*x-1;
        N(:,3) = 0;
        N(:,4) = 0;
        N(:,5) = 0;
        N(:,6) = 4*y;
        N(:,7) = 4*z;
        N(:,8) = -4*y;
        N(:,9) = -4*z;
        N(:,10) = 4*(lambda-x);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1-4*lambda;
        N(:,2) = 0;
        N(:,3) = 4*y-1;
        N(:,4) = 0;
        N(:,5) = 4*z;
        N(:,6) = 4*x;
        N(:,7) = 0;
        N(:,8) = 4*(lambda-y);
        N(:,9) = -4*z;
        N(:,10) = -4*x;
    elseif strcmp(derivative_type, "dz")
        N(:,1) = 1-4*lambda;
        N(:,2) = 0;
        N(:,3) = 0;
        N(:,4) = 4*z-1;
        N(:,5) = 4*y;
        N(:,6) = 0;
        N(:,7) = 4*x;
        N(:,8) = -4*y;
        N(:,9) = 4*(lambda-z);
        N(:,10) = -4*x;
    end
end
end